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Introduction 



Abstract 

We consider a heat conduction model introduced in [2]. This is an open system in 
which particles exchange momentum with a row of (fixed) scatterers. We assume 
simplified bath conditions throughout, and give a qualitative description of the 
dynamics extrapolating from the case of a single particle for which we have a 
fairly clear understanding. The main phenomenon discussed is freezing, or the 
slowing down of particles with time. As particle number is conserved, this means 
fewer collisions per unit time, and less contact with the baths; in other words, the 
conductor becomes less effective. Careful numerical documentation of freezing is 
provided, and a theoretical explanation is proposed. Freezing being an extremely 
slow process, however, the system behaves as though it is in a steady state for 
long durations. Quantities such as energy and fluxes are studied, and are found to 
have curious relationships with particle density. 



1 Introduction 

In this paper, we report on a study of a model introduced in [2] and further studied 
in [3], In this model of 1 -dimensional scattering, there is a row of scatterers, which 
are equally spaced and tied down, the interval between adjacent scatterers being 
called a "gap". Each scatterer carries a single variable, namely its "momentum", 
called q i . Moving about in the system are particles whose momenta and positions 
are denoted by p ■ and x respectively. The particles do not see each other directly, 
and interactions between particles and scatterers follow the rules of elastic scat- 
tering - except that the scatterers do not move. The two ends of the chain are 
connected to heat baths. Unlike [2, 3], which considered the grand canonical case, 
the canonical case is studied in this paper, i.e., the number of particles is fixed, and 
when a particle reaches a bath, a new one is injected immediately with a velocity 
chosen from some distribution. Details of this model are given in Sect. 2.1. 

The fact that the scatterers have no recoil, i.e., their positions are fixed, simpli- 
fies the analysis, making the local dynamical picture more tractable. As we will 
show, locally the dynamics are characterized by rattling: a particle will rattle back 
and forth between two scatterers, carrying momentum from one to the other until a 
certain state is reached, at which time it exits this interval (passing through one of 
the scatterers), and begins to rattle between the two scatterers in its new interval. 

Global problems such as transport and asymptotic states of the chain are beyond 
the reach of rigorous mathematics. Using bath distributions concentrated on two 
short intervals, i.e., injected particles have momenta ~ I L for the left bath and ~ I R 
for the right, we have a fairly complete description of the dynamics in the case of 
a single particle mediating the transfer or momentum among an arbitrary number 
of scatterers. This description is based on a combination of rigorous and heuristic 
arguments, and the results are confirmed in simulations. Extrapolating from this 
single-particle case, we believe we also have a reasonable understanding in the case 
where particle density is low. As particle density increases, the situation becomes 
very complex. We will report on a number of observations and numerical findings, 
some of which we must admit we had not anticipated at the outset. 
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Our original intent was to study the nonequilibrium steady states of this model 
(with Maxwellian bath distributions), which we had thought would resemble steady 
states in similar models such as [5, 8] or [1]; for a review of the broader subject, 
see [7]. What we found instead was that the model considered in this paper freezes, 
i.e., it slows down with time: Through its exchanges with the scatterers, each par- 
ticle acquires, from time to time, very low values of momentum. Once such a 
low momentum is acquired, the particle spends an extremely long time traveling 
between two scatterers, which are unit distance apart. During this time it has no 
influence on the evolution of the system. We find that as time goes on, more and 
more particles are stuck in these low energy states, with fewer and fewer collisions 
occurring per unit time. 1 We further conjecture - not without reason - that "at the 
end of time", the action in the entire system is carried, for the most part, by a single 
particle. That is to say, for arbitrarily long durations, the dynamics are those in the 
single-particle case. 

We will provide numerical documentation of freezing; that takes quite a bit of 
computing time since the process is gradual and very slow. We will also connect 
the no-recoil property of scatterers directly to freezing in one physical dimension, 
thereby providing a theoretical basis for understanding this phenomenon. More 
precisely, we will argue that at least in closed systems (no heat baths), on balance 
particle energy is dissipated through collisions with scatterers that do not recoil. 
In modeling, it is not uncommon to accept an unphysical property to make a model 
more analyzable. In this case we have found that the no-recoil property has an 
unintended consequence. 

While freezing is the main message of this paper, we would argue that quanti- 
ties such as total scatterer energy and fluxes in and out of the chain - observed in 
real time - are entirely relevant. This is because the freezing process is extremely 
slow: We have found that following a relatively brief initial transient, the system 
will settle down to what is perhaps best described as a quasi-stationary state, i.e., 
a very-slowly-varying state which for many purposes can be treated as stationary. 
We have found that during this infinitely long period of quasi-stationarity, the chain 
exhibits macroscopic behaviors that are very challenging to explain. 

This paper is organized as follows: In Section 2, we describe the model and 
explain the mechanics of rattling. Section 3 discusses transport properties of the 
chain assuming quasi-stationarity, and Section 4 is devoted to the documentation 
and discussion of freezing. 

2 Model and local dynamics 

2.1 Model description 

We consider N scatterers at positions % = 1, . . . , N and a "bath" each at positions 
and N + 1. The ith scatterer has a momentum q i G R. There are n particles in 



Freezing here refers to the slowing down of particles, not to falling of scatterer temperatures. On 
the contrary, mean scatterer temperature climbs slowly with time, as will be explained in Sect. 4.3 



4 2.1 Model description 

the system; the momentum of the jth particle is denoted by p £ R. Each particle 
moves with uniform velocity until it reaches either a scatterer or a bath. A collision 
between a particle and a scatterer or bath is called an event in this paper. 

In a particle-scatterer collision, momentum is exchanged as follows: For p = 
Pj and q = q { , the "scattering" replaces these values by p' and q' given by 

p' = -ap + (1 - a)q , 
q' = (1 + a)p + aq . 

In these coordinates, which were used in [2], the scattering matrix is 

—a l — o" 

1 + cr a 



and the energy 



„ 9 1 — <7 9 

E=p 2 + -—q 2 
1 + 0- 



is preserved. In the present paper, it is more convenient to rescale p and q in such 
a way that (p, q) lies on a circle: 



l-o- 



P = P, 9 = 17: — Q 

V 1 + 0" 

In these new coordinates, the matrix takes the form 



-o- Vl^ 



\/r^2 



a" a 



We will be mostly interested in a ~ 1 and thus introduce e = \/l — o- 2 . In these 
coordinates (and we will omit the "hat" from now on), the energy is p 2 + q 2 , and 
the scattering matrix is 



f-y/l - e 2 e 

b ~{ e VT^ 

which geometrically can be seen as a reflection p t-t —p followed by a clockwise 
rotation by an angle 6 = — arctan(e/\/l — e 2 )- 

Once the particle has scattered, it continues on its way, moving to the right of 
the scatterer if p > and to the left if p < 0. 

When a particle reaches a "bath", it disappears and is instantaneously replaced 
by one with a new energy. We have used the word "bath" here to denote an infinite 
source of energy or momenta. No Maxwellian character is assumed, as that will 
be destroyed by the no-recoil property of the scatterers in any case. For simplicity, 
we fix two positive numbers, I L and I R , which are nominal values at the bath, and 
introduce a small amount of randomness by fixing a (small) parameter a > 0. The 
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momenta of injected particles are i.i.d. random variables with a uniform distribu- 
tion on [I L (1— a),I L (l+a)] for the left bath, respectively [— I R (l+a), — I R (1 — a)] 
for the right bath. In general, unless stated otherwise, simulations have been done 
with 

I L = 1 , 4 = 3, e = 0.05, and a = 0.1 , 
and 400 scatterers. It will be convenient to define 

4,eff = 4(! - «) . 4,eff = Jr(1 - «) . 

because, as we shall see, these are the values that occur most prominently when e 
is close to 0. 

2.2 Rattling 

The mechanism by which momentum is transported between adjacent scatterers 
can be summarized as follows: 

(1) Conditions for "crossing". When a particle collides with a scatterer, it can 
be reflected back, or it can pass through the scatterer; the latter is the def- 
inition of "crossing". If a particle with momentum p > collides with a 
scatterer with momentum q < p^f-je, then by the rules of scattering, it will 
be reflected. Similarly for p < and q > p^f-je. In particular, crossing is 
impossible if p ■ q < 0. 

(2) Rattling. Consider an interval bounded by two scatterers with momenta q® 
(on the left) and q® (on the right), and assume that q® > and q® < 0. 
Suppose a particle has just entered this interval. The particle will rattle back 
and forth carrying momentum from one scatterer to the other, and it will exit 
(on one of the two sides) in finite time. After this crossing, the momenta of 
the two scatterers are q 1 w q^ and q 2 ~ (ft. 



2.2.1 Derivation 

We will analyze the case when a particle starts from scatterer 1 , has initial momen- 
tum p > so that it flies to the right, rebounds from scatterer 2, then rebounds from 
scatterer 1 again, and so on. We will follow its trajectory until it exits the interval 
between scatterers 1 and 2, assuming throughout that there is no other particle in 
the picture. 

Applying the matrix S twice (first on p and q 2 and then on p and q ± ) we get the 
sequence: 

V -> ~pV- + e Q2 -> P(V-) 2 - Q2 £ V- + £ Qi » 

Qi ->■ Qi -> qw^-p^^ + s 2 ^ , 
<?2 -+ 12V- + £ p ->• q 2 V- + £ p- 
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300 



Rattling between scatterers and 1 




100 200 300 400 500 

# of events 

Figure 1: Rattling between two consecutive scatterers (for e = 0.01). The horizon- 
tal axis is events. The vertical axis is momentum for g x and q 2 and \p\. (The sign 
of p switches after each collision during rattling.) Note the exchange of the values 
of q 1 and q 2 . 

To discuss this sequence, it is convenient to formulate the problem as a differential 
equation, which takes the form 



£ = ~e(q 2 - qd + 0(e 2 ) , 

an 

^l = +ep + 0{e 2 ). 
an 



(2.1) 



Here n is "discrete time", the passage of each unit of which corresponds to one lap 
(two reflections) for the particle. The solution is then, up to order 2 in e with initial 
conditions p, q t , q 2 : 

p(n) = (q t — q 2 ) sm(env2) / y2 + p cos(£nv2) , 
Qi(n) = \{q x + q 2 ) + \{q x - q 2 ) cos(en\f2) - p s'm(enV2) / V2 , 
q 2 (n) = \{q x + q 2 ) + \{q 2 - q x ) cos(en\/2) - psm(enV2)/V2 . 
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Observe next that if p « then p(n) passes again through for n w ir/(V2£). Thus 
this is about the time when the particle will leave the interval between scatterer 1 
and scatterer 2 (either to the left or to the right). We find from (2.2) and from 
cos(en\/2) = — 1 that qi(n) = q 2 and q 2 (n) = q 1 . Thus, in the case p rs the 
values of q 1 and q 2 are simply exchanged after the rattling, as we have asserted. 

The time evolution of (p(n), (?i(n), q 2 (n)) is illustrated in Fig. 1. 

The question is now on which side the particle will exit. We claim, and have 
checked numerically, that the probability P R to leave to the right as compared to 
that of leaving to the left, P L satisfies approximately 

Pr/Pl = \Ql\/\Qr\ , (2-3) 

where q L is the value of the left-hand q at the time of exit. We have no proof of this, 
but intuitively, this ratio can be understood as follows: For the particle to exit, its 
momentum must satisfy \p\ < e\q\/^/-, with p and q of the same sign. By the study 
above, each collision is a rotation by an angle of order 0(e), in the p, q plane, so 
the probabilities to first satisfy \p\ < e\q L \/-*/- or \p\ < £ \Qr\/\^ are proportional 
to IftJ/lfel- 

2.3 The baths 

We consider next how particles enter the system. For definiteness, consider the 
left bath. Observe that typically, many attempts are made before a new particle 
successfully crosses the first scatterer. This is because in order for a particle to 
leave the chain for (say) the left bath, it must cross the first scatterer from right to 
left. At the time of this crossing, q 1 must necessarily be < 0, and in order for a 
new particle with p w I L to cross it again to enter the chain, q 1 must be raised to 
> I L /e; see the first summary item at the beginning of Sect. 2.2. Raising q 1 to 
this required value is done via a rattling mechanism similar to that discussed in the 
previous subsection, except that bath injections are all « I L . The following gives 
a bound on the number of attempts before a successful entry: 

Lemma 2.1. Suppose (i) particles with momenta p G [p m ; n , P max i are injected, and 
(ii) the first scatterer has initial momentum q and its momentum is not altered by 
any other particle during the period in question. Then after at most 

" (°' ^-Z/2) ) 

collisions with the scatterer the particle will cross into the system, and after cross- 
ing, 

*f < g^ < ^ . (2.5) 

Proof If, before crossing, p > eqf\/\ — e 2 , then the particle will be reflected and 
q changes to q' = q\f\ — e 2 + ep. This can only happen if q < p max \/l — £ 2 /e. 
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Since also p < p max we get 

q' > ep min + q-q(l- \Jl - e 2 ) 



>q + ePmin - PmaxC 1 ~ V^ ~ £ 2 ) Vl ~ &/e 
>q + £(Pmin - Pmax/ 2 ) ■ 



This proves (2.4). Then (2.5) follows at once from q < p max Vi — e 2 /e, the bound 
on p and the scattering rules. □ 

Remark. Particles that get through have p almost exclusively at the low end of 
theoretically admissible range. This is entirely expected: The q's change extremely 
slowly, by 0(s) each hit, while the p's given out by the baths are random lid., so 
the first p to get through is naturally at the low end of the range. 

3 Transport in the chain 

3.1 Single-particle dynamics 

We have a fairly complete description of the dynamics in the case where there is 
a unique particle in the system. As we will see, this simple dynamical system is a 
useful reference point for studying chains with low particle densities (Sect. 3.2). It 
will also be a source of information for a conjectured "final state" for all systems 
(see Sect. 4.4). 

Remark. For simplicity, we will describe all our findings in terms of the numbers 
"1" and "3", which are the injection values, left and right. Our results apply, mutatis 
mutandis, to arbitrary values of I L and I R , though simulations get more difficult 
when the ratio I L /I R is very small or large. 

We find numerically that independently of the initial state of the chain, after 
a long enough transient, all q i become - with small variations - either a "1" or a 
"—3". Here is what we mean: When a particle enters the system, more precisely 
when it crosses its first scatterer, it carries a momentum close to I L eff or — I R eff ; 
see Sect. 2.3. When such a crossing occurs, the scatterer in question generally has 
momentum ~ 1/e times that of the particle; see Sect. 2.2.1. We assert here that 
following a long enough transient, all q i take on essentially these values. When it 
is ~ I L e g/e, we call it a "1"; when it is ~ — I R e g/e, we call it a "—3". This being 
the case, we start the discussion from such a configuration. 

Also, in this single-particle case, the only action is where the particle is, and it 
is simpler to go by events time in lieu of real physical time. This is what is used in 
the discussion to follow. 

Consider the interval between q i and q i+ i- From Sect. 2.2, we see that locally, 
the possibilities are limited: 

(q t , q i+ i) = (1,1): the particle can only enter from the left and exit on the right; 
no change in scatterer energies upon exit; similarly 
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(lij Qi+i) = ( — 3, —3): the particle can only enter from the right and exit on the 
left, with no change in scatterer energies upon exit; 

(lij Qi+i) = (1) — 3): the particle can only enter from the left and exit on the 
right; scatterer energies are transformed from (1, —3) to (—3, 1) 

(<7j, gj + i) = (—3, 1): the particle cannot enter this interval. 

We propose to view all l's as indistinguishable, as are all — 3's, and to consider 
only crossings, i.e., when the particle moves from one interval to another (ignoring 
the rattling that occurs in between). Seen this way, the dynamics of a chain are 
completely described by (i) the walk of the particle, and (ii) flipping or not of 
scatterer energies along the way, where the only "flipping" possible is from (1, —3) 
to (—3, 1). An important observation before going further is that positive momenta, 
i.e., l's, can only move to the right, while negative momenta, i.e., —3's, can only 
move to the left. 

We have found numerically (and have a proof for the single-particle case with 
I L 96 I R ) that a particle cannot stay in an interval forever. We assume for purposes 
of the discussion below that after the particle transforms (1, —3) to (—3, 1) it will 
exit with probability 3/4 to the left and with 1/4 to the right. We now argue that 

in a typical steady-state configuration, there are 3 times as many 1 's 
as there are —3's. 

The argument below is part rigorous and part heuristic, and the phenomenon is 
fully confirmed by numerical experimentation. Notice first the following: 

(a) If we define a "move" to be a swapping of two adjacent energies, then each 
move by a 1 is accompanied by a move by a —3, and vice versa. 

(b) The number of l's that enter the chain is equal to the number of l's that exit (it 
is a steady state) which in turn is equal to the number of —3's that enter (because 
each exit of a 1 is followed by the entrance of a —3) etc. 

In view of (a) and (b) above, we need to show that the —3's move along the 
chain 3 times as fast as the l's. Specifically, we will argue that in each "sequence 
of flips" (to be defined), a —3 moves on average 3 times in a row, while 3 different 
l's move one step each. 

Indeed, consider a configuration of l's and —3's such as 

31111-3111 -31... . 

Suppose the particle is in (1, —3), the middle —3. Then it flips it to (—3, 1), and 
exits the gap. If it goes to the left, it is again in a (1, —3), and the same scenario is 
repeated. This sequence of flips ends when either (i) the particle exits to the right 
after a flip, or (ii) the —3 comes up against the —3 on the left. 

If (i) occurs, we might as well think of the particle as going directly to (1, —3), 
the —3 on the right (because the only exit possible in (1, 1) is to the right) and 
nobody has moved in the meantime. The last paragraph then repeats itself for this 
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^R 


(#q ~ l)/(#g - -/ R ) 


mean p R injected 


exits R / exits L 


3 


2.96 


-2.7304 


1.000031 


5 


5.12 


-4.5508 


1.000019 


7 


7.29 


-6.3712 


1.000006 


9 


9.03 


-8.1915 


1.000005 



Table 1: Illustration of 1 -particle theory. The nominal injections are I L = 1 and 
I R = 3, 5, 7, 9, and a = 0.1. Mean momenta that entered the chain from the right, 



i.e., crossed the rightmost scatterer, are quite close to —I } 



R,eff 



-(1 - a)I R 



— 0.9/ R (third column). The second column confirms our asserted I R : I L ratio of 
I L and — I R in steady states, and the 4th confirms the assertion on left/right exits. 

new —3. If (ii) occurs, then the baton is passed to the —3 on the left and the same 
story happens. 

Now given a configuration of l's and — 3's, imagine doing a walk hopping from 
—3 to —3. In view of our assumption of left/right exits, if the ratio of l's to —3's 
is > 3 : 1, then (i) above will occur more often, resulting in the walker visiting 
the right bath more often than the left. Likewise, when there are too few l's, (ii) 
will happen more often, leading to more frequent visits with the left bath. As noted 
earlier, the two baths are visited with equal frequency in a steady state. Hence in a 
steady state, there are roughly 3 times as many l's as there are —3's. 

Results of numerical verifications for various assertions are shown in Table 1 . 



Summary. In the case of a single particle, we find: 

(1) The system settles down to an events-time steady state in which all scatterer 
energies reflect bath injections with small variations. 

(2) In this steady state, we have 
(i) (qj) = 0, for each i and 

(ii) the energy profile is constant along the chain, with 



(i) 



I/, 



R,effl 



'L,eff 



+ /, 



R,effl 



-Alhn/e 2 ) + 



'L,eff 



-*L,eff + Mr. 



effl 



r(4eff/£ 2 ), (3-D 



(3) Momentum transport along the chain is ballistic in the sense that if we view 
scatterer energies as a sequence of l's and —3's as defined above, then all the l's 
move monotonically from left to right and all the —3's move from right to left. 

3.2 Dynamics of chains with multiple particles 

Consider first the case where there is more than one particle in the chain, but that 
the particle density, g, defined to be the number of particles per scatterer, is <C 1. 
In this case, the dynamical picture can be described as follows: 
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Most of the time, no two particles share a gap or find themselves in adjacent 
gaps, and the actions of the particles are "independent", meaning they do as in 
the single particle case. As a consequence the only prominent changes in scatterer 
energies occur when a particle enters an interval (i, i + 1) with g, > g i+1 ; these 
values are flipped to q[ « q i+1 and q' i+1 « q i as the particle exits this interval. 

When two or more consecutive intervals are occupied by particles, the actions 
of the particle can interfere with one another. One such scenario is illustrated in 
Fig. 2. In the window of time depicted, the first half of the events show particle 1 
flipping q i _ 1 and q^. Before this flip is complete, however, particle 0, which is in 
the interval (i,i + 1), springs into action, causing q i+1 to rise at the expense of the 
pair q i _ 1 and q t , both of which are pushed down to compensate. 

Two particles in two adjacent gaps 
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Figure 2: The interaction of two particles (p (black) and p x (red)) with scatterers 
at locations i — 1, i and i + 1. The horizontal axis is events. The vertical axis is 
momentum for the g, and the \p k \ . Particle rattles between i and i+1 and particle 
1 rattles between i and i — 1. 



While the kind of interaction in Fig. 2 can happen, we have found that in a 
good fraction of encounters, scatterer momenta are in fact swapped in a sequential 
manner, leading to a simple reordering of the g^'s. This can be understood as 
follows: Suppose a particle with momentum p has just entered an interval bounded 
by g 1 w 1/e and q 2 ~ — 3/e. In general, p can be anywhere between and 1 
(roughly speaking), and it takes l/p units of time to traverse the interval. With 
the first reflection, \p\ starts to ramp up by 0(1) per collision, but notice that the 
effects of these initial collisions on the g's are quite insignificant, and that remains 
true until \p\ reaches ~ 1/e, because initially |g| ~ 1/e, and with each collision, it 
changes by only « e\p\ (see Sect. 2.2). The same is true as \p\ ramps down. What 
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this tells us is that in real time, the "flipping" of the q's occurs effectively on a time 
interval that comprises only a fraction of the total "rattling time"; the smaller the 
initial p , the smaller this fraction. Thus even when two particles occupy adjacent 
intervals, if they do not flip the g's at roughly the same time, the action is likely to 
be sequential. 

Baths are much more likely to disrupt this sequential pattern. Recall that it takes 
many attempts before a particle from the bath crosses the first scatterer (Sect. 2.3). 
Consider, for example, the case where particle is in the gap (0, 1) and particle 1 
is in (1, 2), and that particle 1 becomes active while particle is in the middle of 
attempting to gain entrance. Unlike the situation described in the last paragraph, 
the back-and-forth motion of particle between the left bath and scatterer 1 is 
relentless. Its attempt to bring q 1 up can easily interfere with the attempt by particle 
1 to flip q 1 and q 2 , which involves pulling g 1 down. 

Baths have been observed to lead to more complicated scenarios than in Fig. 2. 
We will not enter into such a discussion here, but see Sect. 3.3B. 

The overall picture can be summarized as follows: 

(1) In low particle-density regimes, say for q < 1, one can, to some degree, ex- 
trapolate from the single-particle picture: Many of the q/s are l's and — 3's (see 
Sect. 3.1), and a large majority of the changes in scatterer energies are flips of the 
kind discussed above - except that flips can now occur in different parts of the 
chain in random order. Occasionally, simultaneous action in same or adjacent gaps 
leads to the creation of intermediate g-values. These values, as with all g-values, 
are eventually "flushed out" of the chain. We have found that among the interme- 
diate g-values created, the very small ones (we will call them "zeros") are the most 
stubborn: they tend to remain in the chain for very long times. 2 

(2) As q increases, simultaneous action of particles in neighboring intervals be- 
comes more commonplace. These actions may, in principle, enhance or cancel 
each other, but we see much more of the latter: pulls in opposite directions by 
competing forces lead to decreased amplitudes in the oscillations of q-values. The 
decrease in amplitude on the —3 side is somewhat more pronounced, possibly due 
to the fact that it is easier to interrupt a "longer swing". This phenomenon is im- 
portant; it will help explain some of the observations in the next subsection. 



2 Zeros form barriers which particles cannot cross easily, causing them to oscillate back and forth 
in a short stretch of the chain for a long time. The mechanism is as follows: Assume scatterer 
momenta are 0, 1 or —3, and focus on a particular 0. In order for this to move, a particle must 
enter one of the intervals adjacent to it. If a particle approaches from the left, then the configuration 
is likely to be (1, 0). The particle then flips (1, 0) — > (0, 1) and exits (in all likelihood) to the right, 
our having moved one slot to the left. Similarly, if a particle approaches from the right, then 
(0, —3) —!• (—3, 0), the particle exits to the left, and our moves to the right. Particle approaches can 
occur in any order, but notice that the net move of our cannot exceed the number of particles in the 
chain: Suppose the number of particles is k, and our is k slots to the left of where it was originally. 
The argument above tells us that all k particles must now be to the right of it, so the next approach is 
guaranteed to be from the right! This oscillatory behavior eventually ends when our acquires "by 
accident" a reasonable size (flips do not return exact pre-flip g-values). 
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Figure 3: A snapshot of scatterer momenta q i for 70 scatterers away from the 
boundary of a400-long chain, for 2 different densities g. For g = 0.01, eq i k, 1, — 3 
(injected values) for many i, with a few values near q = 0. For g = 15, eq r values 
are more evenly distributed but concentrated on an interval strictly smaller than 
[—3, 1]; this is due in part to the decreased amplitudes of oscillations explained in 
the text and also to the fact that at any one point in time, more scatterer momenta 
are likely to be in the middle of getting "flipped". 



(3) In higher particle-density regimes, e.g., for q > 5, the process becomes untidy, 
even chaotic. To gain some intuition, we invite the reader to imagine, e.g., for 
g = 10, how 4000 particles move about in a chain with 400 scatterers. These 
particles have quite different velocities, but with a large enough number of them 
present, many will be active simultaneously, and some of the active particles will 
occupy same or adjacent gaps. Each particle behaves as though its mission is to put 
into "the right order" the momenta of the two scatterers at the ends of the interval 
on which it resides, as was discussed earlier. The net effect of all the pulls and tugs 
together will determine the time evolution of the system. 



Snapshots of scatterer energies along a chain for two g's are shown in Fig. 3. 
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3.3 Chain energy and fluxes as functions of g (and time) 

We discuss here two important quantities in relation to energy, namely the mean 
total scatterer energy and influx-outflux rates into the baths. Since fluctuations are 
large, to represent these quantities in a meaningful way, we need to assume there 
is some stationarity - although as we will see in the next section, that in itself is a 
point of contention. Nevertheless, for purposes of the present discussion, we will 
show data from simulations averaged over very large windows of time, and these 
windows will increase in size as time progresses. 

Details on statistics. All of our simulations are performed for a very large number 
of scattering events. To keep measurements comparable for different densities g, 
we define an epoch to be 1.6 • 10 9 • g scattering events. There being 400 scatterers 
in the chains we use, this means that each particle performs on average 4 • 10 6 
scatterings per epoch. Since e = 0.05, rattling between two scatterers usually 
ends after 27r/(n\/2) collisions. This translates into about 200, 000 crossings per 
epoch per particle, which in turn leads us to expect that each particle "sees" the 
bath several hundred times per epoch, the exact number depending on whether 
its motion is closer to ballistic or that of a random walk. Our measurements are 
averaged over each epoch, and we take up to 5000 epochs. 

A. Total scatterer energy. We consider mean total scatterer energy in the chain. 
More precisely, for each value of g we compute the time average of ^ qf in each 
epoch, with data taken over many epochs. The results are shown in Fig. 4. 

Note the agreement with g = (or single-particle) values in Sect. 3.1: substi- 
tuting the simulation values into (3.1) leads to a prediction of (qf) ~ 972, which 
is quite close to the limiting value shown. As g increases, it is evident that total 
scatterer energy decreases. The sharpest decrease occurs for g between and 1. 
We have checked that this is due largely to the accumulation of larger and larger 
numbers of 0's (as explained in Sect. 3.2). As g continues to increase, the de- 
creased amplitudes of oscillations (also explained in Sect. 3.2) will no doubt lead 
to decreased mean scatterer energy, but there is another factor that could potentially 
seriously impact the situation, and that is the rates at which energy flows into and 
out of the chain at the two ends. We now look at these quantities more closely. 

B. Energy fluxes. We consider the following 4 quantities: <J> L in and <£ L out , the 
rates (with respect to physical time) at which energy flows into and out of the 
chain from the left, and <1> R in and <I> R out , the corresponding quantities on the right. 
Notice that once we fix the system parameters N = number of scatterers, g = 
particle density, I L (resp. I R ) = injected momenta from left (resp. right) bath, these 
fluxes are determined entirely by the system. 

Some simulation results are shown in Fig. 5. Observe first that if the system 
is in a steady (or quasi-stationary) state, then by the conservation of energy, one 
would expect <I> L in — $ L out = <!> R out — <1> R in , both quantities being the total flux <!> 
across the system. Fig. 5 does indeed show these two quantities to be equal. The 
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Figure 4: Stalling from a random initial condition, the mean total scatterer energy 
of a chain is measured and plotted as a function of epoch. In time, this energy 
appears to reach some ^-dependent limit, which decreases as g increases. Closer 
inspection, however, shows that except in the case of a single particle, mean energy 
drifts slowly upwards with epoch. 



plot shows that activity per unit time increases with density, which is expected. But 
it also reveals a number of surprising facts. We mention two of them: 

(1) $ as a function of g. The increase is sublinear. To some degree, this can be 
explained by the incomplete flips discussed in Sect. 3.2 and Fig. 3. The inset in 
Fig. 5 shows, in fact, that $ ~ g 1 for 7 « 2/3, which begs for an explanation. 

(2) Contact of system with the right (or hotter) bath. Observe that for g = 10 and 
15, 3> R out/^R i n > 2/3. The system being in a steady state, particles enter and 
leave through the right end of the chain at the same rate, and each entering particle 
carries with it an energy ~ 9 = /| (see Sect. 2.3). This implies that the mean 
energy carried out by exiting particles is > 6, and a particle with p > \/E can exit 
only when q N , the momentum of the rightmost scatterer, is > y/E/e (Sect. 2.2). 
Thus the picture cannot resemble that in Fig. 3. 

We conjecture that the dynamics next to the hotter bath are somewhat volatile: 
Decreased amplitudes in q-oscillations, especially on the negative side, make it 
harder for particles from the right bath to enter the system. As they gather between 
the bath and the iVth scatterer, they must exert a nontrivial downward pull on q N . 
Unable to overcome this pull, a "wall" of large positive q i -values builds up to the 
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Figure 5: The in-flux and out-flux of energy as a function of g. The data show that 
the total flux across the system grows sublinearly with the number of particles. The 
slope of the red line is about 0.68 (inset, log-log). 



left of q N (we have seen this wall many times). The pressures continue to build, 
until at some moment q N swings upwards and a floodgate is opened. 



4 Freezing 



One of the most salient features of this model is that it never seems to approach a 
steady state, in that no matter how long one waits, many quantities continue to drift, 
slightly but perceptibly. We believe this is due to freezing, referring to the slowing 
down of the particles in the system. There is ample mathematical and numerical 
evidence to support this thinking, though we do not have rigorous proofs and it is 
hard for finite-time simulations to provide conclusive evidence for a phenomenon 
that progresses on a slow (logarithmic) scale. In the subsections to follow, we 
examine the problem from several different angles, and attempt to both elucidate 
and document this phenomenon. 
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4.1 Distribution of post-collision particle momenta 

To gain insight into how freezing occurs, we propose to look at distributions of 
particle momenta immediately following collisions with scatterers. We begin with 
a very simple model for which we have a complete description of the dynamics. 

4.1.1 Closed system with one scatterer, one particle, and two walls 

The physical space occupied by this simplified model is the interval [0,2]. At and 
2, there are no baths but two walls, upon contact with which the particle is reflected, 
i.e., p' = —p where p is the momentum of the particle. A single scatterer is placed 
at 1, and the interaction of the particle with the scatterer is as before. Thus the phase 
variables are i] = (q,p, x) where q is the momentum of the scatterer, q 2 + p 2 = c 2 
for some constant c which we may take as 1, and x € [0,2] denotes the position of 
the particle. The flow is denoted by <£ t . 

Let e be the constant in the scattering matrix in Sect. 2.1. 

Theorem 4.1. There is a countable set of e for which <l> t is periodic. These excep- 
tional values of e aside, the following holds: 

(i) For every initial condition rj(0) = (q(0),p(0), x(0)), we have 

^y \ P (t)\dt^o. 

(ii) The expected time between consecutive scattering events is infinite. 

Proof. It is advantageous here to focus on the collision manifold £ = {x = 1}, 
which can be identified with the unit circle {p 2 + q 2 = 1}. We let / : £ — > £ be 
the first return map, and use the following convention: Given 77 = (q,p, 1) G £, 
we first flow, i.e., move right or left according to whether p > or < 0, and do 
the scattering when the particle returns to £. Since S is reflection in p followed by 
rotation by 6 where 6 = — arctan(e/Vl — e 2 ), / is simply rotation by 6. 

As is well known to be the case, if 6 is rational, then /, and hence <J> t , is time 
periodic. This corresponds to a countable set of e. Outside of this exceptional set 
of e-values, / is an irrational rotation, hence it preserves Lebesgue measure m on 
£ and is ergodic. Furthermore, all orbits of / are uniformly distributed on £. 

We now return to the flow $ t , which is what interests us. Let R : £ — > [0, oo) 
be the first return time to S. Then R(q, p) = 2/\p\, so J R dm = oo. It follows that 
for every initial condition 77(0) with no exception, the fraction of time the trajectory 
spends in any given neighborhood of p = tends to 1 as time goes to infinity. □ 

Remark 1. Under the condition that 9 is irrational, the only invariant probability 
measures of this system are {5^ 0x ^,x £ [0, 2]} where 5 denotes point mass at 
the phase point 77. That is to say, all invariant measures are concentrated at phase 
points at which the particle is stationary. This continuum of invariant measures 
taken together can be viewed as a physical measure in the sense that given any 
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initial condition ry(0), ^ J S^dt converges to this family as T — > oo. If one so 
chooses, this can be used as a mathematical definition of freezing. 

Remark 2. The proof of Theorem 4.1 also tells us about the approach to this 
family of stationary measures, i.e., how freezing happens: Fix a very small positive 
number p . The probability of acquiring a value of p with \p\ < p in a collision is 
~ p , which is very small, but once such a value is acquired, the particle becomes 
"inactive" for 2/\p\ units of time. It follows that statistically, the particle spends 
100% of its time in these extremely-low-energy states, punctuated by (rare) periods 
of activity in between. 

4.1.2 Back to the general case 

The analysis of models with N scatterers, n particles and two baths is beyond the 
scope of this paper (or the state of the art of dynamical systems theory for that 
matter). Taking a hint from the simplified model in Sect. 4.1.1, we investigate 
numerically the distribution of momenta following scattering events. 

Ideally, we fix a system, a particle, say particle i, and an arbitrarily chosen ini- 
tial condition. When particle i exits the system, we name the particle that replaces 
it by the same name. We then make a histogram of the momentum of particle i 
after each scattering event that involves it. The nature of the distribution as p — > 
contains much information: If a post-collision momentum is p, then the particle 
will be carrying this value of momentum for the next l/\p\ units of time. Thus 
a PDF that is roughly constant and bounded away from as p — > will imply 
freezing for the same reason as before. 

In practice, it is hard to collect sufficient data for individual particles, so we 
lump all particles together. Results of simulations are shown in Fig. 6. 

Notice that the divergence of the integral J -X dp, which lies at the heart of 
the freezing phenomenon implied here, relies on the one-dimensionality of the 
physical space. On the other hand, had this model been "normal", then when put 
in contact with a Maxwellian heat bath of temperature T, post-collision momenta 
should have a PDF ~ |p|e~^ p where j3 = ^, and should not be ~ 1 as p — > 0. 

4.2 Why this system freezes: a theoretical discussion 

In this subsection we seek to understand why freezing occurs in our model, and 
propose an argument that connects the no-recoil property directly to energy dissi- 
pation in one spatial dimension. Recall that "no recoil" refers to the fact that the 
scatterers do not change their positions (in fact they do not move at all) even though 
they carry "momenta", and the usual momentum exchange between scatterer and 
particle is assumed in collisions. The argument proposed below is heuristic, but we 
believe it sheds light on various features of this model and provides a theoretical 
basis for understanding the freezing phenomenon. 
There are two main ingredients in this argument: 
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Figure 6: Sample distribution of the momenta for epochs 100 and 1800 for 2 densi- 
ties (averaged over snapshots 50 to 100, resp. 900 to 1800). The black line indicates 
the slope corresponding to density l/\p\. Since the time during which a momentum 
p is maintained is l/\p\ this demonstrates numerically that the PDF averaged over 
time would be flat near \p\ = 0. 
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(1) Consider the collision of a particle, whose position and momentum are denoted 
by x and p, with a scatterer, the corresponding coordinates for which are denoted 
by y and q. Had it been normal Newtonian mechanics, the volume form dxdydpdq 
would be preserved, meaning for an infinitesimal volume element V in (x, y,p, q)- 
space, V and & t (V) would have the same volume assuming the collision occurs 
between times and t. Moreover, & t (V) would be stretched in the y-direction 
whenever the scatterer gained energy in a collision. Since assuming the scatterer 
has no recoil is effectively "leaving out" the y-coordinate, it follows intuitively 
that under this assumption, phase volume (in 3D) is contracted when energy is 
transferred from particle to scatterer in a collision. 

This intuition is confirmed in the following direct computation for our model: 

Lemma 4.2. The tangent map through a collision contracts (or expands) the phase 
space volume by a factor \p /p\. 

Proof. Consider a little volume V which is at coordinates (p, q, x) with p < and 
x < initially, and a time t before which all initial conditions in V have made 
exactly one collision with a scatterer located at 0. Then at time t we have 

p' = ~ yfp + m , 

q' =ep + ^q , 

x' =0 + p'(t — \x/p\) = —\f~-p + eq(t — \x/p\) . 

The tangent map is then the matrix 



-s 


e ' 


£ 


V- o 


yf(t ~ \X/ P \) 


—e(t — \x/p\) p' jp 



(4.1) 



and the absolute value of its determinant is \p'/p\. □ 

As we have seen in Sect. 2.2, energy transfer can go either way in collisions, 
so in our model, phase-volume is sometimes expanded and sometimes contracted. 

(2) Next look at the problem from a dynamical systems point of view. For a sys- 
tem that expands phase-volume in some parts of the phase space and contracts it 
in other parts, when Lebesgue measure is transported forward by the dynamics, it 
is likely to accumulate at invariant measures (called physical measures, see [4]) 
that are volume-contracting, or at least nonexpanding, on average. Intuitively, this 
is the reason why trajectories tend to sinks and not sources, attractors rather than 
repellors. It is connected to the idea of entropy production; see e.g., [9]. There are 
also rigorous mathematical results in the same spirit. For example, it is a mathe- 
matical fact that the sum of all Lyapunov exponents of random diffeomorphisms is 
< always, and is = if and only if all the constituent maps preserve the same 
smooth invariant density; see e.g., [6]. 
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We now combine the ideas in (1) and (2). Consider first a closed system, i.e., 
a system that operates in isolation (and is not connected to baths). (2) says that 
asymptotically in time, volume is likely to be decreasing. According to (1), this 
means net energy flow is from particles to scatterers. To summarize, the no-recoil 
property in our model has the following implication: from the view of the parti- 
cles, collisions with scatterers lead — on average — to energy dissipation. This is 
consistent with a slow-down of particles. 

In an open system, however, particles are "recharged" following visits to baths. 
Indeed, the expected number of crossings (referring to the crossing of scatterers) 
a particle makes between visits to the baths is likely to be finite, and the process 
is renewed after each such visit. Thus the freezing process, which occurs for the 
same reason as in closed systems, cannot be completed. This reasoning suggests 
that in open systems the effects of the slow-down should be more severe in longer 
chains. (It also suggests that if the transfer of energy from particles to scatterers 
is strictly positive, then the scatterers will heat up eventually, though we have not 
seen any evidence of that in simulations.) 

4.3 Numerical evidence of freezing 

We provide here two sets of numerical evidence that document the slow-down of 
particles in chains of length 400. 

A. Times between collisions. Perhaps the most direct way to document freezing 
in a chain is to measure the times between collisions as a function of epoch. We 
carried out a study which goes as follows: Start from an initial condition, run the 
system for a very long time, and take snapshots of the system at various points in 
time. At each such time t, we asked how long (in real physical time) it will be 
before X% of the particles have engaged in at least one collision. Suppose this 
happens at time t + T(t, X). Then the slowing down of the system will be reflected 
in the growth of T(t, X) as a function of t for fixed values of X. The results for 
X = 65 for a number of q are shown in Fig. 7. 

B. Energy of a system as function of time. Consider the following two facts: 

(1) Since particles with very low energies do not participate in the evolution of the 
system (other than altering their own positions), the effective density in a system 
should be smaller than its true density q. Moreover, with freezing, one would 
expect this effective density to decrease with time. 

(2) We saw in Sect. 3.3 that the mean total scatterer energy of a system is larger for 
systems with smaller q. 

These two facts together should imply that for a given system with a fixed q, 
its mean total scatterer energy taken over different epochs should creep up slowly. 
Small inclines in the plots in Fig. 4 are indeed evident. 

Important remark. While much of this section is focused on freezing, it is just as 
important to remember that the process occurs extremely slowly: Times between 
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Figure 7: Mean waiting time for 65% of particles as a function of epoch, for several 
densities. Each data point is an average over the epochs from the last data point to 
the current one (and the first from epoch 70 to 100). 
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collisions, as we have seen in Fig. 7, increase only on a log log-scale. Following 
an initial transient, the models considered in this paper remain in a slowly-varying, 
quasi-stationary state for a very long time - indefinitely for practical purposes. For 
this reason, we submit that macroscopic observations such as those in Sect. 3.3 are 
entirely meaningful, even though no well defined limiting values can be reached in 
finite time. 

4.4 The final state 

The following theoretical question begs for an answer, however: What happens as 
time goes to infinity for a chain that is arbitrarily long? In what follows, let us 
consider events time, skipping over the long periods with no collisions occurring 
anywhere in the chain. 

Solo-particle conjecture: For long periods (in events time), the baton is carried 
by a single particle, i.e., it rattles between scatterers and walks about in the chain 
while all other particles are "asleep ". These long solo walks are punctuated by 
(rare) periods during which multiple particles are simultaneously active. The baton 
is passed from particle to particle during these periods of activity. 

The following is the basis of this conjecture: Fix e 1 <C s 2 <C 1. Let us say a 
particle is "inactive" when \p\ < e 1 , "active" when \p\ > s 2 . After the system has 
been running for a long enough time, each particle will be inactive nearly 100% 
of real time, and active for essentially 0% of real time. Unless there are hidden 
correlations (we have no reason to believe there are any), the active periods of a 
particle will most likely occur when all other particles are inactive. Notice also 
that for e 2 small enough, an active particle is likely to remain active for many 
collisions in a row assuming p- values following crossings are more or less random 
(see Sect. 2.2). If the two assumptions above are essentially valid, one will see, in 
events time, long walks by solo particles. 

Observe that during these long walks, the dynamics will, to some degree, re- 
semble those in Sect. 3.1, in that the active particle will be "flipping" scatterer 
energies wherever it goes (though we do not have any basis for speculating if the 
g/s will be mostly — 3's and l's). 

The conjecture above is obviously not testable, nor can any credible evidence 
be produced in reasonable time. We offer some data nevertheless on the distribution 
of lengths of runs, which are defined as follows: Let a system evolve and record the 
sequence of crossings of scatterers as j 1 , j 2 , ■ ■ ■ where j k = i if the kth crossing 
involves particle i. (We have elected to count crossings rather than events, and 
failed attempts to enter the chain are not counted as crossings.) Define a run to be 
a maximal sequence of consecutive j k 's with the same index. Thus a run of length 
n in a system with k particles means that one of the k particles carries out a walk 
that comprises n crossings before a crossing is made by one of the other k — 1 
particles. 

As discussed earlier, we expect these runlengths to get longer with time. A 
numerical experiment with 5 particles was carried out, and some of the results are 



24 



Summary and conclusion 



presented in Fig. 8. The statistics were delicate due to the occurrence of extreme 
events. In the longest runs observed, a single active particle made more than 1.5 • 
10 8 crossings before any of the other 4 particles crossed a scatterer. 




2 3 

log epoch 



Figure 8: Slowing down of particles, for a simulation with 5 particles and 600 
scatterers and the standard heat baths (~ 1, ~ —3). Shown is the percentage of 
crossings spent in runlengths longer than 100, 200, ... as a function of elapsed time 
(on a logarithmic scale, log = 5 corresponding to 8192 epochs of 2 • 10 crossings 
of scatterers, each). The figure is obtained by averaging over 9 realizations.The 
upward trend of the runlength is clearly visible. 



5 Summary and conclusion 

This papers contains an analysis of a model of heat conduction introduced in [2]. 
Our aim was to deduce from the interaction in [2] the large-time dynamical behav- 
ior of the system and to leverage that to shed light on physically relevant quantities 
and phenomena. 

Consider first the simplest case of a single particle. We showed that the particle 
rattles between two adjacent scatterers, systematically moving momentum from 
one to the other - leaving to work on another pair of scatterers after completing its 
task. If one ignores physical time (which is not meaningful here), this translates 
into monotonic or ballistic transport of momenta along chains, accompanied by 
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oscillatory motion in scatterer momenta. 

When more particles are present, the dynamics can be seen as the sum total of 
individual actions, the net effect of a constant tug of war in which each particle 
moves bits of (positive) momentum from the scatterer on its left to the scatterer on 
its right. This way of viewing the microdynamics enables one to explain readily 
certain physical observations that are not obvious otherwise, such as the decrease in 
chain energy as measured by the scatterers as particle density increases. (Reason: 
competing actions of particles in adjacent intervals lead to decreased amplitudes in 
the oscillation of scatterer momenta.) 

But the most striking feature of this model by far is that it freezes, by which 
we refer to the slowing down of the particles. Starting from any initial condition, 
one finds that as time goes on, more and more particles will become inactive, i.e., 
they are in very low energy states, leading to fewer and fewer collisions per unit 
time. On the level of macroscopic observations, this phenomenon manifests itself 
in slowly drifting measurements, independently of the sizes of the epochs used 
in statistical averages. (For example, scatterer energies drift up slowly, consistent 
with the fact that effectively there are fewer particles in the chain.) This drifting, 
however, is very slow, as is the freezing process. We think it is appropriate to view 
the system as being in a quasi-stationary state. 

Finally, we address what we believe is the root cause of the particle slow-down. 
To simplify the local dynamics, the authors of [2] fixed the positions of the scatter- 
ers while retaining the usual rules of energy and momentum exchanges. One way 
to put it is that the scatterers have no recoil. In such a system, phase-space volume 
is not conserved by the dynamics. Specifically, in a particle-scatterer collision, 
phase-space volume is contracted when energy flows from particle to scatterer, ex- 
panded when it flows in the opposite direction. Since dynamical processes have a 
way to head toward volume-contracting regimes (i.e., attractors, and not repellors), 
the tendency here is for particle energy to be dissipated. One is thus reminded 
again that interactions that do not conserve phase volume have consequences. 
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